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We study the shrinking Pearson random walk in two dimensions and greater, in which the direction 
of the A'^*'^ step is random and its length equals A'^^^, with A < 1. As A increases past a critical value 
Ac, the endpoint distribution in two dimensions, P(r), changes from having a global maximum away 
from the origin to being peaked at the origin. The probability distribution for a single coordinate, 
P{x), undergoes a similar transition, but exhibits multiple maxima on a fine length scale for A close to 
Ac- We numerically determine P{r) and P{x) by applying a known algorithm that accurately inverts 
the exact Bessel function product form of the Fourier transform for the probability distributions. 

PACS numbers: 02.50.-r, 05.40.Fb 



I. INTRODUCTION 

In this work, wc investigate the probability distribution of the shrinking Pearson random walk in two and greater 
dimensions, in which the length of the N^^ step equals A^^^, with A < 1. If a walk is at rjv after the A^*'' step, then 
rjv+i is uniformly distributed on the surface of a sphere of radius A^ centered about (Fig.[T]). We assume that the 
walk begins at the origin, and the length of the first step is A" = 1. The random direction for each step corresponds 
to the classic Pearson walk [l], whose solution is well known when the length of each step is fixed. In this case, 
the central limit theorem guarantees that the asymptotic probability distribution of endpoints approaches a Gaussian 
function. 




FIG. 1: Illustration of the first four steps of a shrinking Pearson walk in two dimensions, leading to a displacement r. 



In one dimension, the random walk with exponentially shrinking step lengths exhibits a variety of beautiful prop- 
erties [3, For < 5> the support of the endpoint distribution after N steps, Pn{x), is a Cantor set, while for 
A > ^ the support is the connected interval [— j^j, jz^]- More interestingly, for ^ < A < 1 and for N — )• oo, Pn{x) 
is continuous for almost all values of A, but is fractal on a complementary and infinite discrete set of A values 
A particularly striking special case is A = o = |(\/5 — 1) = 0.618 . . . (the inverse of the golden ratio), where P{x) is 
artistically self-similar on all length scales 0, Q • 

Shrinking random walks in greater than one dimension are much less studied. The probability distribution of 
short Pearson walks with a step size that decays as a power law in the number of steps was treated by Barkai and 
Silbey [9[, while the probability distribution of short Pearson walks with arbitrary unequal step sizes was considered 
by Weiss and Kiefer [l3j . More recently, Rador [ll| studied the moments and various correlations of the probability 
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distribution, and also developed a 1/d expansion method, where d is the spatial dimension, for Pearson walks with 
shrinking steps. 

A physical motivation for this model comes from granular media. If a granular gas is excited and then allowed to 
relax to a static state, the motion of a labeled particle is equivalent to a random walk whose steps lengths decrease 
because of the loss of energy by repeated inelastic collisions. A related example is an inelastic ball that is bouncing on 
a vibrating platform [l2l | , where the velocity of the ball after each bounce essentially experiences a random walk with 
shrinking steps if the vibration is sufficiently weak. Our interest was prompted by M. Bazant [Tsj . who apparently 
introduced the shrinking Pearson walk in an MIT graduate mathematics course on random walks. 

While the distribution of radial displacements, P{r), no longer exhibits self-similar properties, numerical simulations 
indicated that P{r) qualitatively changes shape as a function of A. For A ^ 1, the support of P{r) is confined to 
1 — < r < 1 + Yzj, and the distribution is peaked near r = 1. As A increases beyond ^, the probability of being 
near the origin increases and Pir) eventually exhibits a maximum at the origin when A exceeds a critical value, \c{t). 
For two spatial dimensions, we estimate Ac(r) to be 0.5753882 ± 0.0000003. 

The distribution of a single coordinate, P{x), undergoes a similar shape transition, but at a slight different critical 
value, Ac(a;), that we estimate to be 0.558458 ± 0.000003. More surprisingly, P{x) exhibits up to seven local minima 
and maxima when A ~ Xc{x). The secondary extrema occur on a very fine scale that can be resolved only by a 
high-accuracy numerical method, due to Van Deun and Cools to invert the Fourier transform of the probability 
distribution. 

In the next section, we present some elementary properties of the shrinking Pearson random walk and show how 
to obtain the exact Fourier transform for the radial and single-coordinate probability distributions. In Sec. Ill, we 
apply the Van Deun and Cools algorithm to numerically invert the Fourier transform with high accuracy. From this 
inversion, we outline the behaviors of the radial and single-coordinate probability distributions as a function of A in 
Sec. IV. We briefiy discuss the shrinking Pearson walk in spatial dimensions c? > 2 in Sec. V and conclude in Sec. VI. 



II. BASIC PROPERTIES 



When the length of the A^*^ step decreases exponentially with A^, the shrinking Pearson walk eventually comes 
to a stop at a finite distance from its starting point. Since the direction of successive steps are uncorrclated, the 
mean-square displacement after the A^*'^ step, (r^)^. is given by: 
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In the second line, we use the fact that the directions of different steps are uncorrelated so that the average value of 
all cross terms in the expansion of [ri + r2 + rs + . . . + rjv]^ vanish. We thus obtain the obvious result that (r^) 
grows monotonically with A and diverges as A — > 1, corresponding to the infinite-time limit of the classic Pearson 
random walk. 

Our interest is in the probability distributions of the radial coordinate and a single Cartesian component after A^ 
steps, PNif) and Pn{x), respectively, as well as their A^ — > oo limiting forms, P{r) and P{x). These two distributions 
undergo a transition from being peaked away from the origin for small A, to being peaked at the origin for A greater 
than a critical value. A transition from a unimodal to bimodal probability distribution can be constructed, for 
example, from Brownian motion in media with non- linear shear profiles [Tsj . Here the competition between the flow 
and diffusion drive the transition. In the prsent example, the transition is purely statistical in origin. 

Figure [2] shows the radial distribution for A ~ Ac after a small number of steps to provide a sense for the convergence 
rate to the asymptotic form. For convenience in putting many panels on the same scale, we typically plot the 
distribution P(r)r,„ax versus r/r,„ax, where Tmax = (1 ^ X)~^ is the maximal displacement of the infinite walk. 
Already by A^ = 7 steps, the probability distribution is visually indistinguishable from its asymptotic form. While 
P{r) varies smoothly as a function of A, the position of the global maximum changes dis continuously from being peaked 
at r > to being peaked at r = as A increase beyond a critical value Ac(r). The single-coordinate distribution P{x) 
exhibits a transition from multimodality to unimodality that somewhat resembles the transition for P{r), but is more 
complex in its microscopic details. 

Conventionally, the distribution of the displacment factorizes into a product of single-coordinate distribution, from 
which the radial distribution follows easily. However, in contrast to the classic Pearson walk in which the length of each 
step is the same, the probability distribution for the shrinking Pearson walk no longer factorizes as P(r) — P{x)P{y). 
The differences between the radial and single-coordinate distributions arise because there is a non-trivial correlation 
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FIG. 2: The radial distribution -Pjv(r) for the representative case A — 0.56, with A*" = 1, 2, 3, 4, 5, 6, 7, and 10 steps (upper 
left to lower right). 



between steps in orthogonal directions. If the endpoint of the walk is close to its maximum possible value in, say, the 
x-direction, then the displacement in the y-direction is necessarily small, and vice versa. 

It is worth emphasizing that it is not practical to accurately determine the probability distribution of the Pearson 
random walk with shrinking steps by straightforward simulations. As we shall see, the nature of the transition in 
P{x) is delicate. It would require a prohibitively large number of walks, or a prohibitively fine spatial grid in an exact 
enumeration method, to obtain sufficient accuracy to resolve these subtle features. For this reason, we employ an 
alternative approach that is based on calculating the Fourier transform of the probability distribution — which can 
be done exactly by elementary methods — and then inverting this transform by the highly accurate Van Deun and 
Cools algorithm. 



III. FOURIER TRANSFORM SOLUTION OF THE PROBABILITY DISTRIBUTION 



A. Single-Coordinate Distribution 



We first study the distribution of the (horizontal) x coordinate. To obtain the distribution of x after N steps, 
Pn{x), we start with the Chapman-Kolmogorov equation that relates Pn{x) to Pn-i{x), 

Pn{x) ^ [ dx' PN-i{x')qN{x~x'), (2) 



where qNiw) is the probability of making a displacement whose horizontal component equals w at the N^^^ step. 
Equation ^ states that to reach a point whose horizontal component equals x after N steps, the walk must first 
reach a point with horizontal component x' in iV — 1 steps and then hop from x' to x at the N^^ step. 
We now introduce the Fourier transforms 



P^[k) ^ J dx Pn{x) e'-'"'' , and qM{k) = J dxqN{x)e"'^ , 
to recast the convolution in Eq. ([2]) as the product PN{k) = PN-i{k) qN{k). This equation has the formal solution 

N N 

PN{k)^Po{k)Y[qn{k)^l[qr,{k). (3) 

n=0 n=0 

The latter equality applies for a walk that begins at the origin, so that Paik) = 1. Now qn{x) may be obtained by 
transforming from the uniform distribution of angles to the distribution of the horizontal coordinate in a single step 
by using the relation 

d6 

qn{x) dx = qniO) d9 ^ — , (4) 
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together with x = A" cos 61, to give 
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Although the distribution of angles is uniform, the single-step distribution for the cc-coordinate at the n^^ step has 
a "smile" appearance, with maxima at x = ±A"^^ and a minimum at x = 0. The probability distribution of the 
horizontal coordinate after N steps is a convolution of these smile functions at different spatial scales. It is this 
superposition that gives P{x) its rich properties for A ~ Ac- 

Using the transformation between x and 9 in Eq. ([4]), the Fourier transform of the single-step probability is 

qn{k)^ I dxq^{x)e'''^ l\ee'''^""''°'' ^MkX^-^), (6) 

where Jq is the Bessel function of the first kind of order zero. This result relies on a standard representation of the 
Bessel function as a Fourier integral [l6j . Thus the Fourier transform of the probability distribution in Eq. (jSj may 
be expressed as the finite product of Bessel functions 

^A'(^) = n •^o(^^")- 

To calculate PNix) requires inverting the Fourier transform, 

-1 />oo -[ />oo N—1 

Pn{x) = — dke-'^''PN{k)^ - dfc cosfca- TT Jo(fcA"), (8) 
27r J-oo TT Ja 



where we use the fact that PN{k) is even in k to obtain the second equality. 

Each of the factors Jq in the product in Eq. ([5]) is an oscillatory function of fc, and the product itself oscillates more 
rapidly as the number of terms N increases. The evaluation of integrals with such rapidly oscillating integrands has 
been the subject of considerable research [l9| : in particular, integrals of products of Bessel functions appear in nuclear 



physics [20[, quantum field theory [21|, scattering theory [22|, and speech enhancement software [23|. Recently, Van 



Duen and Cools [14j developed an algorithm that can numerically calculate integrals of power laws multiplied by a 
product of Bessel functions of the first kind quickly and with absolute errors of the order of 10"^^. We use their 
algorithm to compute the probability distribution Pn{x) with this degree of accuracy. To implement their approach, 
we first write [l7| 

AIT / 1 

cosz = -y/— l^— Ji/2{z) - ^/z J3/2(z) 

to express the right-hand side of Eq.® in terms of products of Bessel functions and a power law only. With this 
preliminary, we can directly apply the Van Duen-Cools algorithm to determine PNix) accurately. 



B. Radial Distribution 



For the distribution of the radial coordinate r, P/v(r), we again start with the Chapman-Kolmogorov equation [2| 

Pjv(r) = J dr' Pjv-i(r')Qiv(r - r'), (9) 

where Qn{2) is the probability that the walk makes a vector displacement z at the N^^ step, and we use the angular 
symmetry of the walk to write Pn as a function of only the magnitude of the displacement. Since all angles for the 
iV*'' step are equiprobable, 

Q^(r) = x^^(A^-^-|r|), (10) 
Zirr 

where S{x) is the Dirac delta function. Once again, we use the Fourier transform to reduce the convolution in Eq. (|9]) 
to a product. This recursion has the solution 

N-l N-1 

Pjv(fc) = Po(k) n = n MkX") , (11) 

n=0 n=0 
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with the last equahty appropriate for a walk that starts from the origin. While Eqs. (O and (|TT|) are identical, the 
corresponding distributions in real space are distinct. To obtain Ppf(r), we must calculate 

PN{r) = ^Jdk e-^'^-P^ik) . (12) 

Since P/v (k) is a function of the magnitude of k only, wc can write the integration in polar coordinates and perform 
the angular integration to obtain the spherically symmetric result 

^oo N-l 

Pwir) = — / dkk Mkr) J] Jo(fcA") . (13) 
^ -^^ n=a 

In this Bessel product form, we can again apply the Van Duen-Cools algorithm [li!] to invert this Fourier transform 
numerically. 

IV. THE PROBABILITY DISTRIBUTIONS 

We numerically integrate Eq. (|5|) by the Van Duen-Cools algorithm to give the single-coordinate probability dis- 
tribution Pn{x) whose evolution as a function of A is schematically illustrated in Fig. [3] Notice that there is a value 
A « 0.5567 for which the curvature at the origin vanishes. However, at this value of A the global maximum of the 
P{x) is not at the origin. Thus points where P"{x) = do not help locate the global extrema of the probability 
distribution and we must resort to the numerical integration. 

Since the individual step lengths decay exponentially with N , the finite- A*" distribution PNix) quickly converges to 
its asymptotic — > oo form. For example, for A = 0.56 (close to Xc{x)), the displacement of the walk after 15 steps is 
within 10~^ of its final endpoint. Hence the probability distribution is visually indistinguishable from the asymptotic 
distribution on the scale of the plots in Fig.U] We always use values of N for each A to ensure that xn is within 10~^ 
of its final displacement. For small A, P{x) resembles the smile distribution of the single-step distribution in Eq. ([5]). 



P(X) 




FIG. 3: Schematic and not to scale form of P{x) for increasing A near Ac (bottom to top). For A ~ 0.5567 < Ac the curvature 
at the origin becomes positive, while at A = Ac the location of the maximum in P{x) changes discontinuously. For A — >■ 1, P{x) 
approaches a Gaussian. 

As A approaches Ac from below, the minimum at the origin gradually fills in and disappears for A ~ 0.56. For A > Ac, 
the distribution develops a maximum at the origin that becomes increasingly Gaussian in appearance as A — > 1. 

Unexpectedly, P{x) has multiple tiny maxima near the origin, that are not visible on the scale of Fig. 21 as A 
passes through Ac. The Van Duen-Cools algorithm is essential to obtain sufficient numerical accuracy to observe these 
anomalies. The top line of Fig.[5]shows the quantity Z\(x) = P{x) — 0.387, with the vertical scale magnified by 10'^ to 
expose the minute variations of P{x). At this magnification, one can see the birth of a maximum in P{x) at the origin 
that gradually overtakes the secondary maxima near |a;| ~ 0.2. Consequently, the location of the global maximum of 
P{x) jumps discontinuously from a non-zero value to zero at A = Xc{x) ~ 0.5584558 ± 0.0000003 (as illustrated by 
the middle panel on the top line of Fig. [51 which shows P{x) for Ac — A ~ O (lO~'^) ). 
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FIG. 4: Scaled x-coordinate distribution for the shrinking Pearson walk in two dimensions for A = 0.1, 0.3, 0.4, 0.5, 0.56, and 
0.6 (upper left to lower right). 




FIG. 5: The single-coordinate distribution at highly-magnified scales. Top line: Z\{x) = P{x) ~ 0.387 for A = 0.5565, 0.5584, 
and 0.5598. Bottom line: Ziix) = P{x) - 0.3870562 for A = 0.5564, 0.5567, and 0.5570. 



At a still higher resolution, the nearly flat distribution near x — Q &i magnification 10'^ is actually oscillatory 
at magnification 10^ (Fig. [5] lower line). We see that the small maximum that is born when A passes through 
approximately 0.5565 (Fig. [SJ upper left) actually contains an even smaller dimple that disappears when A > 0.5567 
(middle panel in the lower line of Fig. [5]). To highlight this fine-scale anomaly, we plot, in the lower line of Fig. [SJ the 
quantity Z^ix) = P{x) — 0.3870562 for three A values that are very close to Ac. Intriguingly, we do not find evidence 
of additional anomalous features at a still finer scale of resolution. 

We also use the Van Duen-Cools algorithm to numerically integrate Eq. p3p and determine the radial distribution 
P/v(r). For a small number of steps A^, -P/v(r) changes significantly with each additional step, as was illustrated in 
Fig. [21 Once the number of steps becomes of the order of 10, however, P/v('') is very close to the asymptotic P{r) 
for A ~ Ac(r). The transition behavior in P(r) turns out to be much simpler than that for P{x). For P{r), a peak 
gradually develops at the origin, while the peak r > gradually recedes as A increases. Thus as A passes through 
Ac(r), the location of the global peak of P{r) discontinuously jumps to zero (Fig. We do not find evidence of 
fine-scale anomalies in the radial distribution as A passes through \c{r). 
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FIG. 6: The radial distribution for the shrinking Pearson walk in two dimensions for the cases A = 0.40, 0.50, 0.53, 0.57, 0.59, 
and 0.66. 



V. HIGHER DIMENSIONS 



The approach developed for two dimensions can be straightforwardly extended to higher spatial dimensions. For 
the radial distribution in d dimensions, the single-step distribution QNir) is now 

Q^(-) - ^Jt^J - n ^ (14) 

where fid = 2tt'^/^ /r{d/2) is the surface area of the unit hypersphere in d dimensions and r = |r| is the radial distance. 
The corresponding Fourier transform is [Tsj 

^^^^ ^"d0 sin'^-2 0e''^'^"™^^ ==oi^i(d/2,-fc'A2^/4), (15) 



r(i)r 



^ 2 , 



where o^i(a, z) is the confluent hypergeometric function. The Fourier transform P/v(fc) is then the product of Fourier 
transform of the single-step distributions, and its Fourier inverse gives P/v('r)- By integrating over the d—2 azimuthal 
angles, and then integrating over the polar angle 9, as in Eq. (|15p. the formal solution is 



PN{r) = 1^ J dkk'^-^PNik) £ d9 sin'^-He-'''^'"'', 

2l-d ""o N 



f 

/ dkk''-\F,{d/2,-k\y4) T\oFi{d/2,-k^\^y4). (16) 

"'0 



^'^/2r(d/2) Jo 

Since oFi{i^ + 1, —{z/2)^) oc {z /2)~^ Jy{z) [l^l, we can again numerically determine PNir) by using the Van Duen- 
Cools algorithm. The result of this calculation is that the radial distribution undergoes a second-order transition at 
Ac in which the location of the single maximum continuously decreases to zero as A increases beyond Ac. 

The same formal approach can be used to calculate the distribution P{x). This distribution now remains peaked at 
the origin for all values of A. The physical origin of this property stems from the nature of the single-step distribution. 
The generalization of Eq. ([5]) is 

r , s l{d—3)/2 
g„(x)(X A2("-1)-x2 



This function is flat for d = 3 and peaked at the origin for d > 3. Consequently, the convolution of these single-step 
distributions leads to Pn{x) having a single peak at the origin. 
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VI. DISCUSSION 

We investigated the shrinking Pearson walk, where each step is in a random direction, while the length of the 71*^ 
step is A"~^, with A < 1. Because the step lengths are not identical, one of the defining conditions for the central 
limit theorem is violated. Consequently, there is no reason to expect that the probability distribution for this walk 
is Gaussian. We studied basic properties of the radial probability distribution, P(r), and the distribution of a single 
coordinate, P{x). Because a walk with a large displacement in one direction necessarily implies a small displacement 
in the orthogonal direction, P{r) docs not simply factorize as a product of single-coordinate distributions. The P{r) 
and P{x) are distinct distributions. 

In two dimensions, the radial probability distribution of the shrinking Pearson walk changes from being peaked 
away from the origin to being peaked at the origin as the shrinking factor A increases beyond a critical value Xc{r). 
As this transition in A is passed, the location of the peak changes discontinuously from a non-zero value to r = 0. In 
greater than two dimensions, there is a similar shape transition in the radial distribution, but now the location of the 
only peak goes to zero continuously as A increases beyond Ac(r). 

The single-coordinate distribution P{x) has peculiar features for the specific case of two dimensions. Visually, P{x) 
becomes nearly flat at the origin for A « 0.5565 (middle panel, bottom row of Fig. U). However, at a higher degree 
of magnification, this nearly flat portion of the distribution exhibits fine-scale oscillations, with up to seven local 
extrema. Because additional oscillations can be resolved as the resolution is increased, it is tempting to speculate that 
arbitrarily many oscillations occur at progressively decreasing scales. To test for this possibility, we computed the first 
derivative P}^{x) from Eq. and looked for additional zeros in P'jq{x) as a function of x. Again employing the Van 
Ducn-Cools algorithm, wc find the Pn{x) is strictly positive for x in the range 5 x 10~^ to 10~^ when A = 0.55672, 
but is strictly negative in the same range of x when A = 0.55673. Moreover, P^{x) appears to scale as x^^'^ in the 
range 5 x 10^* < .t < lO^"', so we anticipate no additional zeros for a: — > 0. This numerical test suggests that there 
are no additional oscillations in P{x) beyond those revealed in Fig. \5\ 

We are grateful for financial support from DOE grant DE-FG02-95ER14498 (CAS) and NSF grant DMR0535503 
and DMR0906504 (SR). 
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